En este projecto usaremos el fichero de datos “scrubbed.csv” obtenido en https://www.kaggle.com/datasets/NUFORC/ufo-sightings. Elimininaremos datos de tipo NA y cambiaremos paises vacios (i.e. pais == ““) por”Otros”. Tambien arreglamos el formato de los datos con la libreria lubridate. Usamos la libreria dplyr para eliminar columnas no necesarias.
# Leemos datos de CSV
datos <- read.csv("scrubbed.csv")
# Reemplacamos paises "" por "Otros"
datos$country <- ifelse(datos$country=="","otros",datos$country)
# Formateamos datos
datos$datetime <- as.POSIXct(datos$datetime, format = "%m/%d/%Y %H:%M")
datos$'duration (seconds)' <- as.integer(datos$'duration..seconds.')
datos$latitude <- as.numeric(datos$latitude)
# Eliminamos datos NA
datos <- na.omit(datos)
library(dplyr)
datos <- datos %>% select(-one_of('duration..seconds.', 'duration..hours.min.', 'date.posted'))
AnaliZando los datos por país, lógicamente podemos ver que la mayoría de los OVNIs reportados a la National UFO Reporting Center de EEUU ocurrieron dentro de EEUU. Usaremos la libreria ggplot2 para generar un diagrama representando el numero de OVNIs reportados por país y la librería leaflet para generar un mapa con las longitudes y latitudes.
library(leaflet)
library(ggplot2)
by_country <- aggregate(datos$country, by=list(Country = datos$country), FUN=length)
##Generamos un diagrama representando al número de OVNIs por país
ggplot(by_country, aes(x=by_country$Country, y=by_country$x, fill=by_country$Country)) +
geom_bar(stat="identity", width=1) +
xlab('Paises') +
ylab('Avistamientos') +
guides(fill = guide_legend(title = "Country" ))
leaflet(datos) %>%
addTiles() %>%
addProviderTiles("CartoDB.Positron") %>%
addMarkers(~longitude, ~latitude,
popup = ~state, label = ~city,
clusterOptions = markerClusterOptions())
Por eso, hemos decidido hacer nuestros análisis sobre datos solamente de EEUU. Tomamos una prueba de tamaño 1000.
datos_us <- datos[which(datos$country == 'us'),]
datos_us <- datos_us %>% select(-country)
datos_us <- datos_us[sample(nrow(datos_us), 1000), ]
De la prueba de 1000 datos, se ha hecho el análisis solamente de las variables cuantitativas, las cuales son: (“datetime, latitude, longitude, duration(seconds)”).
Con el comando names(object…) se puede ver todas las columnas del fichero de datos y las referenciadas anteriormente.
names(datos_us)
## [1] "datetime" "city" "state"
## [4] "shape" "comments" "latitude"
## [7] "longitude" "duration (seconds)"
Con summary(object...) podemos ver la media, mediana y
percentiles de varios datos, para calcular la moda usamos la función
mode.
mode <- function(x) {
ux <- unique(x)
ux[which.max(tabulate(match(x, ux)))]
}
summary(datos_us$datetime)
## Min. 1st Qu. Median
## "1946-06-30 21:00:00" "2002-06-20 05:51:30" "2007-07-25 22:00:00"
## Mean 3rd Qu. Max.
## "2005-02-01 17:54:32" "2011-08-20 21:10:30" "2014-05-05 19:30:00"
print(paste("Moda: ", mode(datos_us$datetime)))
## [1] "Moda: 1999-09-01 21:25:00"
summary(datos_us$latitude)
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 21.31 34.21 39.61 38.61 42.27 64.84
print(paste("Moda: ", mode(datos_us$latitude)))
## [1] "Moda: 40.7141667"
summary(datos_us$longitude)
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## -159.37 -116.94 -89.86 -96.00 -81.07 -68.02
print(paste("Moda: ", mode(datos_us$longitude)))
## [1] "Moda: -74.0063889"
summary(datos_us$'duration (seconds)')
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0 45 180 1236 600 259200
print(paste("Moda: ", mode(datos_us$'duration (seconds)')))
## [1] "Moda: 300"
No se ha considerado necesario calcular la varianza y desviación poblacional
| Medida | Formula |
|---|---|
| Cuasivarianza | var(x) |
| Cuasidesviación típica | sd(x) |
| Coeficiente de variación | coef_var(x) |
coef_var <- function(x) {
sd(x) / mean(x)
}
varianza <- var(datos_us$'duration (seconds)')
print(paste("Cuasivarianza: ", varianza))
## [1] "Cuasivarianza: 79582107.6563403"
sd <- sd(datos_us$'duration (seconds)')
print(paste("Cuasidesviación típica: ", sd))
## [1] "Cuasidesviación típica: 8920.88043055955"
cv <- coef_var(datos_us$'duration (seconds)')
print(paste("Coeficiente de variación: ", cv))
## [1] "Coeficiente de variación: 7.21668834470973"
varianza <- var(datos_us$latitude)
print(paste("Cuasivarianza: ", varianza))
## [1] "Cuasivarianza: 30.9979640568702"
sd <- sd(datos_us$latitude)
print(paste("Cuasidesviación típica: ", sd))
## [1] "Cuasidesviación típica: 5.56758152673764"
cv <- coef_var(datos_us$latitude)
print(paste("Coeficiente de variación: ", cv))
## [1] "Coeficiente de variación: 0.144199913362222"
varianza <- var(datos_us$longitude)
print(paste("Cuasivarianza: ", varianza))
## [1] "Cuasivarianza: 346.780377514018"
sd <- sd(datos_us$longitude)
print(paste("Cuasidesviación típica: ", sd))
## [1] "Cuasidesviación típica: 18.6220401007521"
cv <- coef_var(datos_us$longitude)
print(paste("Coeficiente de variación: ", cv))
## [1] "Coeficiente de variación: -0.193986752801377"
Como no se pueden dividir fechas, no calculamos el Coeficiente de variación para datetime.
varianza <- var(datos_us$datetime)
print(paste("Cuasivarianza: ", varianza))
## [1] "Cuasivarianza: 100634912191977392"
sd <- sd(datos_us$datetime)
print(paste("Cuasidesviación típica: ", sd))
## [1] "Cuasidesviación típica: 317230061.929788"
Al calcular la asimetría y curtosis de la latitud determinamos que: Obtenemos una asimetría negativa, por tanto la mayoría de datos se encuentran a la derecha de la media. Obtenemos una curtosis mayor que cero, por tanto, hay una mayor concentración de datos alrededor de la media
library(moments)
x <- (datos_us$latitude)
print(kurtosis(x))
## [1] 3.531213
print(skewness(x))
## [1] -0.170255
print(mean(x))
## [1] 38.61016
hist(x)
### longitude Al calcular la asimetría y curtosis de la latitud
determinamos que: Obtenemos una asimetría negativa, por tanto la
mayoría de datos se encuentran a la derecha de la media. Obtenemos
una curtosis mayor que cero, por tanto, hay una mayor concentración de
datos alrededor de la media. Obtenemos el histograma con los datos de la
longitud.
library(moments)
x <- (datos_us$longitude)
print(paste("Asimetría:", kurtosis(x)))
## [1] "Asimetría: 2.37262590886008"
print(paste("Curtois:", skewness(x)))
## [1] "Curtois: -0.565494647193356"
print(mean(x))
## [1] -95.99645
hist(x)
### datetime Al calcular la asimetría y curtosis de la latitud
determinamos que: Obtenemos una asimetría negativa, por tanto la
mayoría de datos se encuentran a la derecha de la media. Obtenemos
una curtosis mayor que cero, por tanto, hay una mayor concentración de
datos alrededor de la media
library(moments)
x <- (datos_us$datetime)
print(paste("Asimetría:", kurtosis(x)))
## [1] "Asimetría: 11.7306123476029"
print(paste("Curtois:", skewness(x)))
## [1] "Curtois: -2.64878868494963"
hist(x, 50)
### duration (seconds) Al calcular la asimetría y curtosis de la latitud
determinamos que: Obtenemos una asimetría positiva, por tanto la
mayoría de datos se encuentran a la izquierda de la media.
Obtenemos una curtosis mayor que cero, por tanto, hay una mayor
concentración de datos alrededor de la media
library(moments)
x <- (datos_us$'duration (seconds)'
)
print(paste("Asimetría:", kurtosis(x)))
## [1] "Asimetría: 709.306346416546"
print(paste("Curtois:", skewness(x)))
## [1] "Curtois: 25.2127548328487"
print(mean(x))
## [1] 1236.146
Asociados a la Tabla de Frecuencias para algunos datos (excepto histogramas, usados anteriormente para ver mejor las medidas de forma) ### latitud
breaks <- seq(0, 90, by=10)
lat.cumfreq <- c(0,cumsum(table(cut(datos_us$latitude, breaks, right = FALSE))))
plot(breaks, lat.cumfreq,
main="Latitud de los ovnis",
xlab = "Latitud",
ylab = "Latitudes acumuladas")
lines(breaks, lat.cumfreq)
### longitude
breaks <- seq(-200, 100, by=10)
longitude.cumfreq <- c(0,cumsum(table(cut(datos_us$longitude, breaks, right = FALSE))))
plot(breaks, longitude.cumfreq,
main="Longitud de los ovnis",
xlab = "Longitud",
ylab = "Longitudes acumuladas")
lines(breaks, longitude.cumfreq)
Si quieres tratar de avistar un UFO, la siguiente información que te proporcionamos te será de suma utilidad, como el país, estado y ciudad donde debes ir para realizar un avistamiento, además de la forma en la que debes buscar y en las horas a las que debes buscarlo,también tenemos otros datos curiosos como la duración y los años que fueron más probables de avistar.
print("Probabilidad de avistar un ovni por país:")
## [1] "Probabilidad de avistar un ovni por país:"
dat <- read.csv("scrubbed.csv")
dat_country<-aggregate(dat$country, by=list(Country = dat$country), FUN=length)
dat$country <- ifelse(dat$country=="","otros",dat$country)
numAvistamientos<-sum(dat_country$x)
proobau<-dat_country[which(dat_country$Country == 'au'),]
probbauval<-(proobau$x/numAvistamientos)*100
print(paste("Australia: ",probbauval,"%"))
## [1] "Australia: 0.669720659264054 %"
proobau<-dat_country[which(dat_country$Country == 'ca'),]
probbauval<-(proobau$x/numAvistamientos)*100
print(paste("Canadá: ",probbauval,"%"))
## [1] "Canadá: 3.73450181745755 %"
proobau<-dat_country[which(dat_country$Country == 'de'),]
probbauval<-(proobau$x/numAvistamientos)*100
print(paste("Alemania: ",probbauval,"%"))
## [1] "Alemania: 0.130707563611014 %"
proobau<-dat_country[which(dat_country$Country == 'gb'),]
probbauval<-(proobau$x/numAvistamientos)*100
print(paste("Gran Bretaña: ",probbauval,"%"))
## [1] "Gran Bretaña: 2.37140865408555 %"
proobau<-dat_country[which(dat_country$Country == 'us'),]
probbauval<-(proobau$x/numAvistamientos)*100
print(paste("Estados Unidos: ",probbauval,"%"))
## [1] "Estados Unidos: 81.056117113977 %"
proobau<-dat_country[which(dat_country$Country == ''),]
probbauval<-(proobau$x/numAvistamientos)*100
print(paste("otros países: ",probbauval,"%"))
## [1] "otros países: 12.0375441916048 %"
maxcountry<-sum(by_country$x)
ggplot(by_country, aes(x=Country, y=(x/maxcountry)*100, fill=Country)) +
geom_bar(stat="identity", width=1) +
xlab('Paises') +
ylab('Probabilidad') +
guides(fill = guide_legend(title = "Country" ))
Debido a que Estados Unidos es donde hay mayor probabilidad de avistar un UFO, vamos a calcular los estados donde es mas probable lograrlo
datos_us$state <- ifelse(datos_us$state=="","otros",datos_us$state)
by_state <- aggregate(datos_us$state, by=list(State = datos_us$state), FUN=length)
avistTot<-sum(by_state$x)
mostseekstate<-by_state[which(by_state$x>16),]
maxstate<-by_state[which(by_state$x == max(by_state$x)),]
print(paste("El estado con mayor probabilidad de avistar un ovni es",maxstate$State,"con una probabilidad de un ",(maxstate$x/avistTot)*100,"%"))
## [1] "El estado con mayor probabilidad de avistar un ovni es ca con una probabilidad de un 13 %"
print(paste("A continuación se muestra una gráfica con los estados donde mayor probabilidad hay de avistar un ovni"))
## [1] "A continuación se muestra una gráfica con los estados donde mayor probabilidad hay de avistar un ovni"
ggplot(mostseekstate, aes(x=mostseekstate$State, y=(mostseekstate$x/avistTot)*100, fill=mostseekstate$State)) +
geom_bar(stat="identity", width=1) +
xlab('') +
ylab('Probabilidad') +
theme(axis.text.x=element_blank())+
guides(fill = guide_legend(title = "Estados" ))
Ya que hemos calculado los estados, haremos lo mismo con las
ciudades
by_city <- aggregate(datos_us$city, by=list(City = datos_us$city), FUN=length)
avistTot<-sum(by_city$x)
mostseek<-by_city[which(by_city$x>4),]
maxcity<-by_city[which(by_city$x == max(by_city$x)),]
print(paste("La ciudad con mayor probabilidad de avistar un ovni es",maxcity$City,"con una probabilidad de un ",(maxcity$x/avistTot)*100,"%"))
## [1] "La ciudad con mayor probabilidad de avistar un ovni es chicago con una probabilidad de un 0.8 %"
## [2] "La ciudad con mayor probabilidad de avistar un ovni es las vegas con una probabilidad de un 0.8 %"
## [3] "La ciudad con mayor probabilidad de avistar un ovni es seattle con una probabilidad de un 0.8 %"
print(paste("A continuación se muestra una gráfica con las ciudades donde mayor probabilidad hay de avistar un ovni"))
## [1] "A continuación se muestra una gráfica con las ciudades donde mayor probabilidad hay de avistar un ovni"
ggplot(mostseek, aes(x=City, y=(x/avistTot)*100, fill=City)) +
geom_bar(stat="identity", width=1) +
xlab('') +
ylab('Probabilidad') +
theme(axis.text.x=element_blank())+
guides(fill = guide_legend(title = "Ciudades" ))
En cuanto a la forma que tienes que buscar en el cielo para lograr un avistamiento exitoso
by_shape <- aggregate(datos_us$shape, by=list(shape = datos_us$shape), FUN=length)
by_shape$shape <- ifelse(by_shape$shape=="","otros",by_shape$shape)
max_shape<-sum(by_shape$x)
maxshape<-by_shape[which(by_shape$x == max(by_shape$x)),]
minshape<-by_shape[which(by_shape$x == min(by_shape$x)),]
print(paste("El ovni mas probabilidad de ser avistado son aquellos ovnis con forma",maxshape$shape,"con una probabilidad del",(maxshape$x/avistTot)*100,"%"))
## [1] "El ovni mas probabilidad de ser avistado son aquellos ovnis con forma light con una probabilidad del 22.3 %"
print(paste("Mientras que los más raros con mayor dificultad de avistamiento son aquellos con forma"))
## [1] "Mientras que los más raros con mayor dificultad de avistamiento son aquellos con forma"
print (paste(minshape$shape))
## [1] "cross"
print(paste("con una probabilidad de avistamiento del",(minshape$x[1]/avistTot)*100,"%"))
## [1] "con una probabilidad de avistamiento del 0.3 %"
ggplot(by_shape, aes(x=by_shape$shape,y=(by_shape$x/max_shape)*100,fill=by_shape$shape)) +
geom_bar(stat="identity", width=1) +
xlab('') +
ylab('Avistamientos') +
theme(axis.text.x=element_blank())+
guides(fill = guide_legend(title = "Formas" ))
La hora también es importante, hay que saber cuando buscar y cuando
descansar
by_horas<- aggregate(datos_us$datetime, by=list(Hour = format(datos_us$datetime,"%H")), FUN=length)
maxhora<-by_horas[which(by_horas$x == max(by_horas$x)),]
maxhour<-sum(by_horas$x)
print(paste("La hora en la que es más probable avistar un ovni son las",maxhora$Hour,':00',"con una probabilidad del",(maxhora$x/avistTot)*100,"%"))
## [1] "La hora en la que es más probable avistar un ovni son las 21 :00 con una probabilidad del 16.1 %"
plot(by_horas$Hour,(by_horas$x/maxhour)*100,type="s",col="dark red", xlab = "Horas del día", ylab = "Probabilidad", main = "Probabilidad de avistamiento por hora")
También es interesante ver cuando hubo una mayor probabilidad de avistar un ovni
anios<- aggregate(datos_us$datetime, by=list(Date = format(datos_us$datetime,"%Y")), FUN=length)
maxanio<-anios[which(anios$x == max(anios$x)),]
maxyear<-sum(anios$x)
print(paste("El año con mayor probabilidad en caso de avistamiento fué",maxanio$Date,"con una probabilidad del",(maxanio$x/avistTot)*100,"%"))
## [1] "El año con mayor probabilidad en caso de avistamiento fué 2013 con una probabilidad del 9.3 %"
plot(anios$Date,(anios$x/maxyear)*100,type="l",col="dark blue", xlab = "Linea temporal", ylab = "Probabilidad", main = "Probabilidad de avistamientos por año")
Observamos la probabilidad de lograr un avistamiento de larga
duración
durability<- aggregate(datos_us$`duration (seconds)`, by=list(duracion = datos_us$`duration (seconds)`), FUN=length)
maxduracion<-sum(durability$x)
maxdur<-durability[which(durability$x>1),]
durationProb<-durability[which(durability$x<=10),]
sumaprobduracion<-sum(durationProb$x)
print(paste("La probabilidad de ver un UFO durante menos de 10 segundos es del",(sumaprobduracion/maxduracion)*100,"%"))
## [1] "La probabilidad de ver un UFO durante menos de 10 segundos es del 11 %"
durationProb<-durability[which(durability$x<=60),]
sumaprobduracion<-sum(durationProb$x)
print(paste("La probabilidad de ver un UFO durante menos de un minuto es del",(sumaprobduracion/maxduracion)*100,"%"))
## [1] "La probabilidad de ver un UFO durante menos de un minuto es del 56.5 %"
durationProb<-durability[which(durability$x>60),]
sumaprobduracion<-sum(durationProb$x)
print(paste("La probabilidad de ver un UFO durante mas de un minuto es del",(sumaprobduracion/maxduracion)*100,"%"))
## [1] "La probabilidad de ver un UFO durante mas de un minuto es del 43.5 %"
plot(maxdur$duracion,(maxdur$x/sumaprobduracion)*100,type="l",col="blue", xlab = "Segundos", ylab = "Probabilidad", main = "Probabilidad de avistamientos por más de un segundo")
A continuación realizaremos algunas operaciones con sucesos de
probabilidad
prob1 <- datos_us[which(datos_us$city == 'seattle'&format(datos_us$datetime,"%H")>12),]
prob2<-aggregate(prob1$city,by=list(city=prob1$city),FUN=length)
print(paste("La probabilidad de ver un UFO en seattle por la tarde",(sum(prob2$x)/avistTot)*100,"%"))
## [1] "La probabilidad de ver un UFO en seattle por la tarde 0.7 %"
prob1 <- datos_us[which(datos_us$city == 'seattle'|format(datos_us$datetime,"%H")==21),]
prob2<-aggregate(prob1$city,by=list(city=prob1$city),FUN=length)
print(paste("La probabilidad de ver un UFO en seattle o verlo a las 21:00 es del",(sum(prob2$x)/avistTot)*100,"%"))
## [1] "La probabilidad de ver un UFO en seattle o verlo a las 21:00 es del 16.8 %"
prob1 <- datos_us[which(datos_us$city != 'seattle'&format(datos_us$datetime,"%H")!=21),]
prob2<-aggregate(prob1$city,by=list(city=prob1$city),FUN=length)
print(paste("La probabilidad de ver un UFO en una ciudad que no sea seattle a una hora distinta a las 21:00 es del",(sum(prob2$x)/avistTot)*100,"%"))
## [1] "La probabilidad de ver un UFO en una ciudad que no sea seattle a una hora distinta a las 21:00 es del 83.2 %"
prob1 <- datos_us[which(datos_us$city != 'seattle'&format(datos_us$datetime,"%H")==21),]
prob2<-aggregate(prob1$city,by=list(city=prob1$city),FUN=length)
print(paste("La probabilidad de ver un UFO en una ciudad que no sea seattle a las 21:00 es del",(sum(prob2$x)/avistTot)*100,"%"))
## [1] "La probabilidad de ver un UFO en una ciudad que no sea seattle a las 21:00 es del 16 %"
prob1 <- datos_us[which(datos_us$city != 'seattle'&format(datos_us$datetime,"%H")==21&(datos_us$shape == 'light'|datos_us$shape=='unknown')),]
prob2<-aggregate(prob1$city,by=list(city=prob1$city),FUN=length)
print(paste("La probabilidad de ver un UFO en una ciudad distinta de seattle a las 21:00 con una forma que sea lumínica o desconocida",(sum(prob2$x)/avistTot)*100,"%"))
## [1] "La probabilidad de ver un UFO en una ciudad distinta de seattle a las 21:00 con una forma que sea lumínica o desconocida 5 %"
prob1 <- datos_us[which(datos_us$city == 'seattle'&format(datos_us$datetime,"%H")!=21&(datos_us$shape != 'light'&datos_us$shape!='unknown')&datos_us$`duration (seconds)`>10),]
prob2<-aggregate(prob1$city,by=list(city=prob1$city),FUN=length)
print(paste("La probabilidad de ver un UFO en seattle a una hora distinta de las 21:00 con una forma que sea distinta de lumínica o desconocida pero con una duracion de mas de 10 segundos",(sum(prob2$x)/avistTot)*100,"%"))
## [1] "La probabilidad de ver un UFO en seattle a una hora distinta de las 21:00 con una forma que sea distinta de lumínica o desconocida pero con una duracion de mas de 10 segundos 0.5 %"
Respondemos algunas preguntas de probabilidad condicional
maxdat<-count(datos)
prob1 <- datos_us[which(format(datos_us$datetime,"%H")==15),]
prob2<-count(prob1)
print(paste("Dado que son las 15:00 que probabilidad hay de ver un UFO?",(prob2/avistTot)*100,"%"))
## [1] "Dado que son las 15:00 que probabilidad hay de ver un UFO? 1.6 %"
prob1 <- datos_us[which(datos_us$city == 'richmond'),]
prob2<-count(prob1)
print(paste("Si vivo en Richmond que probabilidad tengo de ver un UFO?",(prob2/avistTot)*100,"%"))
## [1] "Si vivo en Richmond que probabilidad tengo de ver un UFO? 0.1 %"
prob1 <- datos_us[which(datos_us$shape=='changing'),]
prob2<-count(prob1)
print(paste("En caso de ver un ovni que probabilidad hay de que sea cambiante?",(prob2/avistTot)*100,"%"))
## [1] "En caso de ver un ovni que probabilidad hay de que sea cambiante? 2.8 %"
prob1 <- datos[which(datos$country != 'us'),]
prob2<-count(prob1)
print(paste("Si vivo fuera de Estados Unidos que probabilidad hay de avistar un ovni?",(prob2/maxdat)*100,"%"))
## [1] "Si vivo fuera de Estados Unidos que probabilidad hay de avistar un ovni? 18.9425722359854 %"
Como los datos, excepto los extremos, están aproximadamente sobre la linea, podemos decir que exluyendo los extremos, los datos siguen una distribucion normal.
dia <-as.numeric(format(datos_us$datetime,'%d'))
qqnorm(dia, pch = 1, frame = FALSE)
qqline(dia, col = "steelblue", lwd = 2)
También tenemos la media, moda, varianza y desviación típica de estos
datos.
print(paste('Media:',mean(dia)))
## [1] "Media: 15.033"
print(paste('Moda:',mode(dia)))
## [1] "Moda: 1"
print(paste('Varianza:',var(dia)))
## [1] "Varianza: 76.8667777777778"
print(paste('Desviación típica:',sd(dia)))
## [1] "Desviación típica: 8.76737006050148"
Función de ditribución acumulada
f1 <- pnorm(dia,mean(dia),sd(dia))
plot(dia,f1,main = "Funcion de distribucion acumulada")
Calculamos la probabilidad de que se avisten ovnis durante los primeros 10 dias de cada mes de todos los años.
P(X≤10)
print(paste('La probabilidad es:',100 * pnorm(10,mean(dia),sd(dia)),"%"))
## [1] "La probabilidad es: 28.2963463300533 %"
La probabilidad de que se avisten ovnis durante los ultimos 5 dias de cada mes
P(X>25)
print(paste('La probabilidad es:',100 * (1-pnorm(25,mean(dia),sd(dia))),"%"))
## [1] "La probabilidad es: 12.7804901649419 %"
La probabilidad de que se avisten ovnis entre los dias 10 y 20 de cada mes.
P(10≤X≤20)
print(paste('La probabilidad es:',100 * (pnorm(20,mean(dia),sd(dia)) - pnorm(10,mean(dia),sd(dia))),"%"))
## [1] "La probabilidad es: 43.1520610555124 %"
Podemos representar este último gráficamente:
regionX=seq(10,20,0.01)
xP <- c(10,regionX,20)
yP <- c(0,dnorm(regionX,15,12),0)
curve(dnorm(x,15,12),xlim=c(0,30),yaxs="i",ylim=c(0,0.035),ylab="f(x)",
main='Densidad P(10<X<20)')
polygon(xP,yP,col="orange1")
box()
Podemos averiguar cuantos ovnis fueron divisados en 2012. Para ello, decimos nuestro landa es: “número de ovnis avistados en 2012”.
year <-as.numeric(format(datos_us$datetime,'%y'))
landa <- sum(year == 12)
plot(dpois(0:200, landa), type = "h", lwd = 2,
main = "Función de masa de probabilidad",
ylab = "P(X = x)", xlab = "Número de ovnis avistados en 2012")
Una vez calculada la tasa de llegada de los ovnis en 2012 la distribucion de Poisson, podemos dibujar la gráfica de la función de densidad exponencial.
# Rejilla del eje X
x <- seq(0, 0.1, 0.01)
# lambda
plot(x, dexp(x, landa), type = "l",
ylab = "f(x)", lwd = 2, col = "red")
De la misma forma tenemos la funcion de distribución de exponencial
acumulada
# Rejilla de valores del eje X
x <- seq(0,0.1, 0.01)
pexp(x,landa)
## [1] 0.0000000 0.5934303 0.8347011 0.9327945 0.9726763 0.9888910 0.9954834
## [8] 0.9981637 0.9992534 0.9996965 0.9998766
# lambda
plot(x, pexp(x, landa), type = "l",
ylab = "F(x)", lwd = 2, col = "orange")
La probabilidad de que se divisen menos de 100 ovnis en 2012 P(X<=100) es:
ppois(100,landa)
## [1] 0.8650998
Probabilidad de que se divisen entre 80 y 95 ovnis en 2012 P(80< X< 95)
A modo ilustrativo, podemos representarlo en una gráfica
ppois(95,landa) - ppois(80,landa)
## [1] 0.5646824
pois_sum <- function(lambda, lb, ub, col = 4, lwd = 1, ...) {
x <- 0:(lambda + lambda * 2)
if (missing(lb)) {
lb <- min(x)
}
if (missing(ub)) {
ub <- max(x)
}
plot(dpois(x, lambda = lambda), type = "h", lwd = lwd, ...)
if(lb == min(x) & ub == max(x)) {
color <- col
} else {
color <- rep(1, length(x))
color[(lb + 1):ub ] <- col
}
lines(dpois(x, lambda = lambda), type = "h",
col = color, lwd = lwd, ...)
}
pois_sum(landa,80, 95, lwd = 2,
col = 2, ylab = "P(X = x)", xlab = "Ovnis avistados en 2012")
Gráfico de la función de distribución exponencial
plot(ppois(40:100, landa), type = "s", lwd = 2,
main = "Función de distribución exponencial",
xlab = "Número de eventos", ylab = "F(x)")
También podemos obtener los cuantiles correspondientes de la distribucion exponencial
plot(qpois(seq(0,1,0.001), landa),
main = "Función cuantil",
ylab = "Q(p)", xlab = "Cuantiles",
type = "s", col = 3, xaxt = "n")
axis(1, labels = seq(0, 1, 0.1), at = 0:10 * 100)
Para los datos, vamos a calcular la probabilidad de se vean ovnis el dia 24 de cada mes.
Con lo que si se ven ovnis los días 24, son aciertos, y al contrario son fallos.
dia <-as.numeric(format(datos_us$datetime,'%d'))
tam<- 1000
x <- (dia == 24)
y <- sum(x)
probEnsayo <- (y/tam)
print(probEnsayo)
## [1] 0.025
Funcion de probabilidad binomial
plot(dbinom(1:100, tam, probEnsayo), type = "h", lwd = 2,
main = "Función de probabilidad binomial",
ylab = "P(X = x)", xlab = "Número de éxitos")
Probabilidad de encontrar a lo largo de 69 años (2014 - 1945) mas de 35
ovnis el día 24 de cada mes:
P(X > 35)
1 - pbinom(35, tam, probEnsayo)
## [1] 0.02104234
Probabilidad de encontrar a lo largo de 69 años (2014 - 1945) menos de 25 ovnis el día 24 de cada mes:
P(X <= 25)
pbinom(25, tam, probEnsayo)
## [1] 0.5529257
Ahora graficmente la probabilidadd de encontrar entre 25 y 35 ovnis.
P(25 <= X <= 35)
binom_sum <- function(size, prob, lb, ub, col = 4, lwd = 1, ...) {
x <- 0:size
if (missing(lb)) {
lb <- min(x)
}
if (missing(ub)) {
ub <- max(x)
}
plot(dbinom(x, size = size, prob = prob), type = "h", lwd = lwd, ...)
if(lb == min(x) & ub == max(x)) {
color <- col
} else {
color <- rep(1, length(x))
color[(lb + 1):ub ] <- col
}
lines(dbinom(x, size = size, prob = prob), type = "h",
col = color, lwd = lwd, ...)
}
binom_sum(0 :100, probEnsayo, lb = 25, ub = 35, lwd = 2,
ylab = "P(X = x)", xlab = "Número de éxitos")
Funcion de distribucion binomial
plot(pbinom(10:50, tam, probEnsayo), type = "s", lwd = 2,
main = "Función de distribución binomial",
xlab = "Número de éxitos", ylab = "F(x)")
Hemos usado la distribucion binomial para el teorema central del límite,
Realizamos el histograma de las 6 últimas muestras y de la media muestral,
tamMuestra=1000
numMuestras=25
#hemo
#Se deja ese hueco vacío aunque salga warning
mat=matrix(, nrow = numMuestras, ncol = tamMuestra)
mediax=vector()
# Genera numMuestras muestras aleatorias
# Para cada muestra, almacena la media en el vector anterior
for (i in 1:numMuestras){
x=rbinom(tamMuestra,200,probEnsayo)
mat[i,]=x
mediax[i]=mean(x)
}
# Muestra la media de las 6 últimas muestras generadas
print(c(mediax[numMuestras],mediax[numMuestras-1],mediax[numMuestras-2],
mediax[numMuestras-3],mediax[numMuestras-4],mediax[numMuestras-5]))
## [1] 5.166 5.037 5.058 5.025 5.042 5.007
# Muestra el histograma de las 6 últimas muestras y de la media muestral
par(mfrow=c(3, 3))
hist(mat[numMuestras,],probability=TRUE)
hist(mat[numMuestras-1,],probability=TRUE)
hist(mat[numMuestras-2,],probability=TRUE)
hist(mat[numMuestras-3,],probability=TRUE)
hist(mat[numMuestras-4,],probability=TRUE)
hist(mat[numMuestras-5,],probability=TRUE)
hist(mediax,probability=TRUE)
Vamos a construir un intervalo de confianza del 95% para la duracion media de los avistamientos en EEUU.
#filtramos por duracion menor o igual a 7200, para descartar valores muy altos que son inútiles en nuestro análisis
retval <- subset(datos_us, datos_us$`duration (seconds)` < 7201, select = c(state, `duration (seconds)`))
intervalo <- t.test(retval$'duration (seconds)', conf.level = 0.95)#sacamos el intervalo de confianza
print(intervalo$conf.int)#intervalo de confianza
## [1] 599.1089 766.4317
## attr(,"conf.level")
## [1] 0.95
print(mean(retval$'duration (seconds)'))#media
## [1] 682.7703
print(paste("El resultado indica que la media de la variable duracion (seconds) es de ", mean(retval$duration..seconds.), " el cual se encuentra con una confianza del 95% en el intervalo ", intervalo$conf.int[1], " ", intervalo$conf.int[2]))
## [1] "El resultado indica que la media de la variable duracion (seconds) es de NA el cual se encuentra con una confianza del 95% en el intervalo 599.108912262483 766.431738144021"
car::qqPlot(retval$'duration (seconds)', pch=19,
main='QQplot para la duración de avistamientos',
xlab='Cuantiles teóricos',
ylab='Cuantiles muestrales')
## [1] 5 87
hist(retval$'duration (seconds)', freq=TRUE,
main='Histograma para la duración de avistamientos',
xlab='Duracion (s)',
ylab='Frecuencia')
#Podemos observar que la mayoría de avistamientos duran muy poco tiempo.
Dado que la mayoría de los datos son de EEUU, queremos saber si los OVNIs realmente prefieren estar ahí. Para esto, haremos un contraste de la duración del avistamiento entre EEUU y todos los demás. Podemos asumir que como los datos son del National UFO Reporting Center (NUFORC) basado EEUU, datos de otros países no son datos poblacionales.
Basaremos el estudio en las siguentes hipótesis:
| Hipótesis | Contraste |
|---|---|
| H0 (Nula) | µ(EEUU) = µ(Otros) |
| H1 (Alternativa) | µ(EEUU) != µ(Otros) |
Donde µ(EEUU) y µ(Otros) son la duración media de un avistamiento en EEUU y en otros países respectivamente. Tomaremos el nivel de significancia α = 0.1
alpha <- 0.1
datos_us <- datos[which(datos$country == 'us'),]
datos_us <- datos_us %>% select(-country)
datos_otros <- datos[which(datos$country != 'us'),]
datos_otros <- datos_otros %>% select(-country)
# Eliminamos valores extremos
duracion_us <- datos_us$`duration (seconds)`[!(datos_us$`duration (seconds)` %in% boxplot.stats(datos_us$`duration (seconds)`)$out)]
duracion_otros <- datos_otros$`duration (seconds)`[!(datos_otros$`duration (seconds)` %in% boxplot.stats(datos_otros$`duration (seconds)`)$out)]
Como tenemos un número de datos muy grande, usamos el teorema central del límite (verificando que es distribución normal) para calcular µ(EEUU). Dado que las duraciones son muy variables, asumimos que las desviaciones desconocidas son distintas.
muestra_duracion_us <- replicate(100, sample(duracion_us, size = 1000, replace = TRUE) %>% mean())
car::qqPlot(muestra_duracion_us)
## [1] 20 48
var_us <- var(muestra_duracion_us)
x_us <- mean(muestra_duracion_us)
Hacemos lo mismo para las duraciones de otros países
muestra_duracion_otros <- replicate(100, sample(duracion_otros, size = 1000, replace = TRUE) %>% mean())
car::qqPlot(muestra_duracion_otros)
## [1] 100 5
var_otros <- var(muestra_duracion_otros)
x_otros <- mean(muestra_duracion_otros)
Calculamos la region de rechazo
f <- ((var_us / 100 + var_otros / 100) ^ 2) / ((var_us/100)^2/99 + (var_otros/100)^2/99)
f <- round(f)
cuantil <- qt(c(.025, .975), df=f)
print(paste("La region de rechazo es T < ", cuantil[1], " o T > ", cuantil[2]))
## [1] "La region de rechazo es T < -1.9720790337785 o T > 1.9720790337785"
estadistico <- abs(x_us - x_otros) / sqrt(var_us / 100 + var_otros / 100)
print(estadistico)
## [1] 6.013512
Como nuestro estadistico está en la region de rechazo, rechazamos la hipótesis nula. Eso significa que efectivamente, la duración de un avistamiento es mayor en EEUU.
Podemos verificar esto usando el test t de student ya que nuestra n > 30
t.test(x=muestra_duracion_us, y=muestra_duracion_otros, alternative = "two.sided", conf.level = 0.9)
##
## Welch Two Sample t-test
##
## data: muestra_duracion_us and muestra_duracion_otros
## t = 6.0135, df = 197.42, p-value = 8.649e-09
## alternative hypothesis: true difference in means is not equal to 0
## 90 percent confidence interval:
## 6.29813 11.07161
## sample estimates:
## mean of x mean of y
## 266.0416 257.3567